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ABSTRACT 

The finite element methods (FEM) have proved to be a powerful technique 
for the solution of boundary value problems associated with partial 
differential equations of either elliptic, parabolic, or hyperbolic type. 
They also have a good potential for utilization on parallel computers 
particularly in relation to the concept of domain decomposition. 

This report is intended as an introduction to the FEM for the 
nonspecialist. It contains a survey which is totally nonexhaustive, and it 
also contains as an illustration, a report on some new results concerning two 
specific applications, namely a free boundary fluid-structure interaction 
problem and the Euler equations for inviscid flows. 


This work was supported under the National Aeronautics and Space 
Administration under NASA Contract No. NAS1-18107 while the author was in 
residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), Hampton, VA 23665-5225. 


i 



INTRODUCTION 


It is totally impossible to survey the theory of finite element methods 
within a few pages, and the object of this article is to describe for the 
nonspecialist some very basic ideas and concepts in finite elements approxima- 
tions and to discuss some future trends in the theory without any attempt at 
being exhaustive. Beside this survey part, this article contains in Sections 
8 and 9 a report on some new results concerning two specific problems, viz. ji 
free boundary fluid structure interaction problem and the Euler equations for 
inviscid flows. 

There is no agreement about the first appearance of the method. Finite 
element methods have probably been used for many years for computing and 
engineering purposes in a more or less explicit form. R. Courant mentions in 
[10] the approximation of a function in IF? by continuous piecewise linear 
functions on a triangulation, and this may be the first appearance in the 
mathematical literature. 

Although it is difficult to track the first appearance of the method, 
there is no doubt that the first systematic and large scale utilizations of 
the finite element methods (FEM) occurred in the sixties in solid mechanics 
engineering. The period coincides of course with the first computers and the 
early stages of what we now call scientific computing. Probably the reasons 
which made FEM immediately popular among solid mechanics engineers is that, as 
we recall below, the foundations of the FEM coincide with some very funda- 
mental concepts in solid mechanics. The method has then spread with different 
levels of response in fluid mechanics, optimization and control theory, and 


among mathematicians 
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Alike the solid mechanists, mathematicians (numerical analysts and some 
more theoretically oriented mathematicians) have been working in FEM because 
the methods are appropriate for mathematical treatment and are very close in 
their fundamental concepts to the ideas and tools which are used in the 
mathematical treatment of the linear and nonlinear boundary value problems by 
functional analysis. 

The mathematical and engineering literature on FEM for partial 
differential equations is abundant, and there is no way to survey it here. 
The questions that we address are the following ones: In Section 1, we recall 
the principle of weak formulations, and in Section 2, we recall the role of 
domain decomposition in the context of structural mechanics. We return to 
domain decomposition in Section 9 as it relates to future developments in the 
FEM in relation with parallel computation and some possible extensions of the 
method. Sections 3 to 5 are devoted to the description of very typical 
mathematical results. Section 3 describes the general mathematical framework 
and the most common finite elements. Section 4 provides some convergence and 
error results while Section 5 is an introduction to mixed and hybrid finite 
elements. Some specific applications (among many others) of the FEM are then 
described. Section 6 is related to the Navier-Stokes equations. Section 7 
deals with fluid-structure interactions problems, and Section 8 deals with the 
applications of FEM to the solution of the Euler equations. Finally, as 
indicated above, we return in Section 9 to domain decomposition and the role 
that this can play in future developments for FEM. 
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FOUNDATIONS OF THE FEM 

The finite element methods lie on two fundamental ideas: 

- the weak formulation of a boundary value problem, 

- the domain decomposition , i.e., the decomposition of the domain 
corresponding to the problem into smaller subdomains, the elements. 

As mentioned above, both ideas are closely related to basic concepts of 
solid mechanics. The weak formulation of a boundary value problem coincides 
with the virtual work theorems and energy principles in the statics of 
solids. Domain decomposition is also an extrapolation of the natural approach 
in structural mechanics where large structures consist of smaller substruc- 
tures which are properly connected or assembled, and the study of the large 
structure is reduced to that of the elementary structures and their 
connections. 

1. Weak Formulations 

We begin by recalling briefly the weak formulation of some boundary value 
problems in solid and fluid mechanics. Other examples of weak formulations 
will appear below (abstract boundary value problems). 

l.a. Weak Formulations in Solid Mechanics 

Consider a solid body which fills at rest a region Q of IF? with 
boundary T. We assume that the body is subjected to volumic forces of 

density f = (^^ 2 ^ 3 ) in Q and to surface (traction) forces of surface 

density F = ( F 1> F 2 ,F 3^ on sorae part F 1 of r anci reac ^ es a new 

equilibrium position. The unknowns of the problem are : 

- the field displacements, u = (u^^,^), u(x) x£Q, representing the 
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displacemeat between the position at rest of a particle xefi and its 
new equilibrium position x + u(x). 

- the boundary stress tensor field, a = (a^). 

Under the assumption of small displacements, the equilibrium equations 

read 


( 1 . 1 ) 


( 1 . 2 ) 


3 3o. . 

I 

j=l j 


3 l '» ^ 


in ft 


on fj. 


where v = (v^.v^.v^) is the unit outward normal on r. 

Usually the displacement u is given on the complementary part 


r 0 of r,.r 0 -rv 1 


(1.3) 


u = U on Tq. 


The so-called set of statically admissible stress tensors S a ^(f,F) is 

the set of tensor fields a satisfying (1.1) and (1.2). The set C (U) 

ad 

is the set of kinematically admissible displacement fields, i.e., the set of 
u's satisfying (1.3). 

The equations (1.1) - (1.3) which hold for any material are supplemented 
by the constitutive equations of the material which depend on the material and 
connect stresses and displacements. Without describing these relations, we 
can already see the weak formulation of the problem. Let o, u be solution 
of (1.1) - (1.3) and let v be another kinematically admissible field of dis- 
placements, v eC ad (U) (and w = v - u€C ad = C a< j^ 0 ^‘ We (1*1) by 
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w.p add these relations for i = 1,2,3, integrate over Q, and use Green 
formula and (1.2) (1.3). We obtain 


(1.4) 


f o e (w)dx = / f iWi dx + / F^dr, for all w eC 1 


ad’ 


where the Einstein summation convention has been used and e(w) 
is the strain tensor 


(e^ (w)) 


. 3w. 3w. 

e ^ (w) = 2 ( 3TT + dir h 


If we remember that a = a(u) due to the constitutive law, we find 
that (1.4) is the weak formulation for the displacements 1 . For instance, in 
the simplest case of linear elasticity we have pointwise 


(1.5) 


a^(x) = A i j kl e kl (u)(x), for all xea, 


where the coefficients ^ijkl define a linear positive invertible operator 
A in the space of symmetric tensors of order two. Whence (1.4) becomes 

(1.6)/ A i j kl (u)e kl (u)e 1 ^(w)dx =/ f^dx + / f^dr, for all weC^. 

$ r ^ 

In linear and nonlinear elasticity, the weak formulation (1.4) (or (1.6)) 
coincides with the relation given by the virtual work theorem. It leads also 
to energy principles. 


1 


A similar formulation is available for the stresses a, 


l 
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l.b. Weak Formulations in Fluid Mechanics 

Weak formulations in fluid mechanics do not have a physical interpreta- 
tion as natural as in solid mechanics. They have been introduced by J. Leray 
([16], [17], [18]) for the study of weak (i.e., nonregular) solutions of the 

Navier-Stokes equations in an attempt to explain turbulence by the appearance 
of singularities in the curl vector of the flow. Although we do not know yet 
if such singularities arise in space dimension three, there is no doubt that 
the contribution of J. Leray has been a fundamental step for the mathematical 
treatment of the Navier-Stokes equations by the modern methods of functional 
analysis and also for the numerical treatment of the equations in Computa- 
tional Fluid Dynamics. 

Consider for example the Navier-Stokes equations of an incompressible 
fluid in the stationary case. The fluid fills a bounded region ft of 1^ 
with boundary T. In the Eulerian representation of the flow, the unknowns 
are the velocity field u = (uj,^,^) and the pressure field p; u = u(x), 
x€ft is the velocity of the particle of fluid at x, and p(x) is the 
pressure at point x. We have the equations 

(1.7) -vAu + (u*V)u 4* Vp = f in ft 

(1.8) div u « 0 in ft , 

where v > 0 is the kinematic viscosity and f represents volumic 

forces. Equation (1.7) is the equation of conservation of momentum. Equation 

(1.8) is the incompressibility equation, i.e., the equation of mass conserva- 
tion. If f is materialized and moving with velocity U, then the nonslip 



condition on r is 


(1.9) 


u = U on r. 


Let l/(U) be the space of functions satisfying (1.8) and (1.9). 
u£l/(U) and if v is a test function in (/(U) , w = v - u£ 1/(0). We 
scalar product of (1.7) with w (pointwise in IE?), integrate over 
use Green's formula. We have 

3 u i 3w^ 

-/ Au wdx = -/ Au iWi dx = / — — dx 

ft ft £2 J J 

/ grad p»wdx = / pw»vdr -/ p div wdx = 0. 

ft r £2 

Whence p disappears and we obtain the weak formulation of (1.7) - (1.9): 


Then 
take the 
ft, and 


( 1 . 10 ) 



and for every w£((0) 


3u. 9 wj 3 

l sir 57: dx + . I l 

ft j j i,J = l ft 


3u i 

u i nr "j dx 

i J 


3 

= 1 / 

i=l 


fi w idx. 


It is equivalent to say that u satisfies (1.10) or that u satisfies (1.7) 
- (1.9). The striking fact in formulation (1.10) is that the pressure 

disappears and we are left with an equation involving u only. Once u is 
found we know from mathematical results that there exists p which is defined 
up to an additive constant by (1.7). However, in the practice of numerical 
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computations p is obtained differently, in general, as the Lagrange multi- 
plier of the constraint div u = 0 (see Section 5 below). 

2. Domain Decomposition 

In structural mechanics it is natural to compute a complicated structure 
by considering the smaller substructures of which it is made. Each substruc- 
ture is well modeled, its behavior is well understood, and then the mechanical 
engineers model the interaction (contact laws, etc.) of the different 
components to obtain the description of the full structure. 

As mentioned before, finite elements in solid mechanics have started as 
an extrapolation of this idea to continuous bodies: the full solid body is 
decomposed into smaller elements; a simplified constitutive law is adopted on 
each element; and a simplified version of the constitutive law leads to 
simplified interactions laws between the contiguous elements. 



Figure 2.1 
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Similarly, the particle and cells methods in fluid mechanics which are 
very close to the finite element methods are based on a simplified analysis of 
the flow in small cells with simplified fluid transfer laws. The generaliza- 
tion and mathematization of the FEM have led to a more systematic view and a 
more systematic approach. 

Beside discretization, there are several other good reasons to decompose 
a large domain into smaller subdomains. These reasons are also at the heart 
of future developments in scientific computation and probably finite 
elements. We will return on this important question in Section 8. 


MAIN METHODS - MAIN MATHEMATICAL RESULTS 

We give an overview of some typical finite element methods and some 
typical mathematical results which have been obtained. 

3. The Usual Finite Elements Methods 


3. a. A Model Problem 

We consider as a model problem the following mathematical problem. 

We denote by 12 an open bounded domain of I?, with boundary T, and 
we consider a Laplace equation, 

(3.1) -Au + u = f in 12, 


with associated boundary conditions of Dirichlet and Neuman type 
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(3.2) 


u = 0 on Tq 


(3.3) 


du 

3v 


0 on r x 


where is a partition of f. In the two-limit cases 

r Q = r, r L = 0 and r Q = 0, r x = r we obtain respectively the Dirichlet and 

Neuman problems; the general case is a mixed boundary value problem. 

Let V be the space of functions u satisfying (3.2) and possessing a 
certain level of regularity which we do not specify at the moment. The solu- 
tion u of (3.1) - (3.3) belongs to V and if v is a test function in V, 
we multiply (3.1) by v, integrate over ft, and apply Green's formula. 

Thanks to (3.3) (and v = 0 on Tq) we find 


(3.4) 


n 

-J Au vdx = £ 
ft i=l 


f fLH 
ft 9x i 



dx 


and thus 


(3.5) 


u£V and 

a(u,v) = (f,v), for all v£V 


where 


(3.6) 


a(u,v) = I f dx + / uvdx 

i=l ft i i ft 


and 
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(3.7) (f,v) = / f(x) v(x)dx 

n 

Is the scalar product in 

Conversely, it can be proved (under suitable regularity assumptions) that 
if u satisfies (3.5) then u is the solution of (3.1) - (3.3). Equation 
(3.1) is derived from (3.5) by appropriate methods using distribution deriva- 
tives; (3.2) follows from "u€V" whereas (3.3) is a boundary condition hidden 
in (3.5). This is a general fact with weak formulations like (3.5): some 
boundary conditions of the problem are contained in the definition of the 
space V, and some boundary conditions are contained in the equation (3.5). 

Let us give a more precise definition of the space V. Roughly speaking, 
the space V will be the space of all functions u vanishing on Tq and 
such that a(u,u) < ». More precisely it is easy to see that the expression 

{ a(u,u)} ' l 

is a norm on the space of continuously differentiable functions on TT which 
vanish on r^. We define V as the completion of this space for this norm; 
we obtain the space 

(3.8) V = {veH^n), v| = 0} 

0 

where H*(ft) is the Sobolev space 

H^a) = {veL 2 (n), !^-eL 2 (n), 


(3.9) 


i" 1 y • • > n } 



- 12 - 


More generally H m (f2), the Sobolev space of order m, is the space of 

2 

functions u square integrable in J2 (u€L (£2)) such that all derivatives 
of order m are square integrable also. 

3.b. Abstract Boundary Value Problem 

The situation in (3.5) is typical of many linear elliptic boundary value 
problems. The abstract setting is the following one: 

- We are given a Hilbert space V (norm II • II y) and a bilinear form a 
on V x V which is continuous, i.e., 

(3.10) There exists M < » such that 

a(u,v) < M II ull Itvlly, for all u,v€V 

and coercive, i.e., 

(3.11) There exists a > 0 such that 

a(u,u) a Ilully2, for all u€V. 

- We are given also a linear continuous form Jt on V, i.e., an element 
of the dual V' of V 

and then the problem is 


(3.12) 


To find u€V such that 
a(u,v) = <Jl,v>, for all v€V 
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Despite its simplicity, (3.12) is applicable to many interesting boundary 
value problems in mechanics and physics. The existence and uniqueness of a 
solution u of (3.12) is classically provided by the Lax-Milgram Theorem (see 
for instance [28]). 

More generally nonlinear elliptic boundary value problems can be set in a 
form similar to (3.12) if we allow V to be a Banach space and a to be 
nonlinear with respect to its first argument (i.e., a maps V x V into Rand 
is linear with respect to its second argument). For instance, it follows 
readily from (1 . 10) that the stationary Navier-Stokes equations (1.7) - (1.9) 
with U = 0 can be written in this form. Similarly, consider the problem 
(3.1) - (3.3) and replace the linear equation 

-Au + u = f in SI 


by the nonlinear one 

(3.13) -Au + p(u) = f in SI, 

where p is a polynomial of odd degree with a positive leading coefficient. 
Then (3.13), (3.2), and (3.3) can be set in a form similar to (3.12) 

V = {v€H 1 (n)nL a+1 (fl), v | _ = 0}, 

0 

a(u,v) =lf dx + / p(u)vdx, 

i=l SI i i SI 

where o is the degree of the polynomial p (see [28]). 
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In the nonlinear case, there are no general assumptions on a covering 
all the interesting situations, and we will restrict ourselves to specific 
examples . 

3,c« General Form of Finite Element Approximations 

The discretization of the abstract boundary value problem (3.12) consists 
in choosing 

- a family ^ V h^h€H fin * te dimensional approximations of V 

- a family (a^u^, v h^h€H of bilinear forms on x V h which 
approximate a. 

Roughly speaking there are two types of discretizations produced by the 
finite elements: 

- the conforming finite elements in which the are subspaces of V of 

higher and higher dimensions as the parameter h* 0, 

- the nonconforming finite elements in which the are not subspaces 

of V. 

Of course finite elements have been only used in space dimension n = 2 
and at a less developed stage when n = 3. We consider first the case where 
ft is a polygonal set. A basic ingredient of FEM is a triangulation of 
ft. By this, we mean a suitable covering of ft by either 

- a family of triangles, 

- a family of rectangles whose sides are parallel to the axes (or more 
general quadrilateral sets), 

- or a combination of triangles and rectangles (or quadrilateral sets). 

The triangles or rectangles are the (finite) “elements. 11 The space 


consists of functions of a given type (usually a polynomial) on each element 
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which are properly connected. The values of the functions of or their 

derivatives at some particular points of the elements (vertices, mid- 
edges,...) are the nodal values which fully determine the functions in V^. A 
natural basis of consists of the so-called shape functions ; These are 

the functions of whose nodal values are 1 for one of them and 0 for 

all the others. In most cases these functions have a ,, sraall ,, support, and 
this leads to fairly sparse matrices for the discretized problem. When a 
function v is defined on ft (or on an element K) , its interpolant on 
ft (or K) denoted r^v (or r^v) is the function of (or the elementary 

function on K) which has the same nodal values as v. 

3.d. Conforming Finite Elements (n = 2) . 

For second order elliptic boundary value problems, the basic space V 
is H*(ft) or a product of such spaces or a subspace of such spaces. 

The simplest and most common elements used in this case are the 
elements on triangles and the Q| elements on rectangles. ?i (respectively 
P n ) is the set of polynomials of degree 1 (resp. ^ra), whereas Qi (resp. 
Q m ) is the set of polynomials of degree 1 (resp. _< m) with respect to each 
variable. 

Some other typical elements used for second order boundary value problems 
are depicted in Figure 3.1. We will return to the and elements 

after we briefly describe the elements in Figure 3.1. 

Triangles 

linear: polynomials of degree <_ 1 on the triangles, nodal values - 


values at vertices 



Linear, quadratic, cubic triangles 



Linear, quadratic, cubic rectangles 



Conforming Reduced Cubic Reduced quadratic 

Herrnite elements Triangle Triangle 

Figure 3.1: Conforming Finite Elements (n = 2) 
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quadratic: polynomials of degree <,2 on the triangles; nodal values = 

values at vertices and midedges. 


cubic: polynomials of degree <^3 on triangles; nodal values = 

values at vertices, barycenter and 1/3 points on edges. 


reduced cubic: polynomials of degree 3, vanishing at the barycenter on 

each triangle; nodal values = values at vertices and 1/3 
points on edges. 

Rectangles 

linear: polynomials of degree £ 1 in each variable on 

rectangles; nodal values = values at vertices. 

quadratic: polynomials of degree <^2 in each variable; nodal 

values = values at vertices, midedges, and center. 


cubic: polynomials of degree £ 3 in each variable; nodal 

values = values of function at 16 different points (see 
Figure 3.1). 


reduced quadratic: polynomials of degree 2 in each variable satisfying a 

linear relation (on each rectangle); nodal values « 
values of function at vertices and midedges. 
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All functions obtained by these elements are globally C® (continuous) 
except the quadratic Hermite triangle which produces C 1 approximants 
(continuously differentiable functions). 

More special elements can be found in the literature; cf. for instance 
the book of P. G. Ciarlet [9] on the mathematical side, and the book of 
Zienkiewicz [33], the work of Argyris [1] and others on the engineering or 
mechanical sides. 

The more sophisticated elements produce better (more precise) results but 
need more computing time and a good expertise in finite elements technology. 
In a nonspecialized industrial environment, the tendency seems to be the 
utilization of simple elements of degree one or at most two with a suitable 
refinement of the mesh. 

As mentioned above the simplest and most commonly used elements are the 

elements on triangles and the elements on rectangles with sides 

parallel to the x and y axes. Let us mention also the quadrilateral 
elements described hereafter. 

a 

Let K denote the square (0,1) x (0,1) in the 5, n plane. We 
observe that a mapping F with Qj components 

! a + b£ + cq + d£ q 
a' + b'£ + c'q + d'£q 

A 

can map K on any arbitrary quadrilateral K of the x,y plane. The 

image by F of a line in the £ ,n plane is generally a curved line of the 

x,y plane. However the lines x = constant, y = constant, and in particular 

* 

the boundary of K are mapped by F onto straight lines of the x,y 
plane. A natural element on the quadrilateral K is the image by F~* of 
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the Qj 

I 




Figure 3.2 

element on K: 

q = q(5,n)eQ 1 qoF 1 = qoF ^x.y). 

In general these elements are not polynomials on K. They are, however, 
easy to use and their explicit expression is rarely used. 

3.e. Conforming Finite Elements (n = 3) 

The triangulation is now the covering of ft (- a polygonal set) by 
either tetrahedrons or 3-D rectangles whose edges are parallel to the axes or 
combinations of those. 

The most common elements are the 

- linear, quadratic, and cubic tetrahedrons. 

- linear, quadratic, and cubic 3-D rectangles. 

The definitions of these elements are the same as above in the two-dimensional 
case replacing triangle by tetrahedron and rectangle by 3-D rectangle. For 




Linear, quadratic, cubic tetrahedrons 


Linear, quadratic, cubic rectangles 


Figure 3.3: Conforming Finite Elements (n 
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the cubic (tetrahedron and 3-D rectangle) elements, the nodal values are shown 
on Figure 3.3, All these elements lead to functions which are globally C° 
(continuous) but not more. 

3.f. Nonconforming Finite Elements 

As indicated above, nonconforming finite elements produce approximate 
function spaces which are not subspaces of V. For instance, the 

linear nonconforming triangle described below produces, when applied to 
Problem (3.12), approximate functions which are highly discontinuous. Still, 
it may be useful to use such elements in at least two cases: 

- Fluid flow problems where, due to the incompressibility condition divu 
= 0, the linear triangle elements cannot be used in a straightforward 
manner. 

- Higher order problems, like the biharmonic problem where most elements 

described above fail to produce C* functions, and thus the approximate 

2 

spaces would not be included in H (ft) (= the natural space for a 

biharraonic problem) . 

Nonconforming linear elements 


2-D Case (Triangles) 

Polynomials of degree 1 on the triangles 
Nodal values = values at midedges 
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3-D Case (Tetrahedron) 

Polynoms of degree <^1 on the tetrahedrons 
Nodal values « values at barycenter of faces 

The global functions are totally discontinuous with discontinuities along 
the edges of triangles (or faces of tetrahedrons) except for the barycenters 
(of edges or faces). The method is nevertheless convergent and efficient, 
particularly for fluid flows: see the book of F. Thomasset [32] which is 

fully devoted to the utilization of these elements in 2-D flows. 

3. g. Curved boundaries 

Curved boundaries can be approximated by polygonal lines. Alternatively 
one can use the so-called isoparametric elements: the element is the image by 

an appropriate (simple) mapping of a triangle or a rectangle and the function 
reduces on the element to the composition of a polynomial with that mapping. 
A similar situation occurred with the quadrilateral elements. 

4. Convergence and Error Estimate 

Concerning convergence and error estimates the situation is different for 
linear and nonlinear problems. 

4. a. Linear Problems 

Two type of results have been derived in relation with error computation 
and convergence (see for instance P. G. Ciarlet [9]): 

- interpolation error, 


- approximation error 
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When is a conforming finite element space and u is a function in V (or 

usually in a smaller space), we consider the interpolant r^u of u in 
(this is the finite element function which assumes the same nodal values as 
u, whereas r^u is the interpolant of u on an element K); the interpola- 
tion results give an upper bound of the norm of u - r^u in V and other 
spaces. The approximation results are of a different nature: when u£V is a 
solution of a problem like (3.5) and u^€V^ is a so ^ ut ^ on t ^ ie associated 
discrete problem, then the error between u and u^ is estimated for various 
norms. In the optimal cases the error between u and u^ is of the same 
order as that of the distance of u to V^. 

The general results are too abstract to be presented in detail here; we 
will just recapitulate the error estimates corresponding to the elements 
described above. 

For an element K let p v denote the radius of the smallest ball 

containing K, let p' denote the radius of the largest ball included in 

K. 

K, and let = P K /p'. 

The analysis is made under the assumptions that 


(4.1) 


p h 


= Sup p„ -*■ 
' K£T h 


0 


and 

o u = Sup a remains bounded from above. 

K eT u K 

h 

is a function in V and r^V its interpolated function in V^, 
the H m semi-norm of v - r^v on an element K of the triangu- 


(4.2) 

If v 

we consider 
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lation T h and on the whole domain ft : 

i v - vL.k - j. ] l I D “ (V - v)i 2d 4 l/2 

V la J =m K ) 

iv ■ v.* ■ Lu i idV ' v,l S 1/2 ' 


where D a is a partial derivative of order [a ] = m and the sum is 
extended to all such derivatives. 

For the elements described above, the interpolation result is the 
following one: 

On an element K assume that the interpolation operator r K is such 
that r p = p for each polynom p of degree k, and assume that 
r Is linear continuous from H^ + 1(K) into H m (K), 0 < m < k + 1. Then 



We can also assemble the results on the different elements K 
angulation and obtain a similar bound on all of Q (when 

polygon fully covered by the elements): 


I \ y n k+l-m mi i 

l v • vim, a i c p h °h |vl lc+l,(! - 


k+1 

for all veH (a), 


of a tri- 
U is a 
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Finally in dimensions 2 or 3 

- for the linear elements (triangles, rectangles, tetrahedrons, 3-D 

rectangles) if: 

k = 1, v€H 2 (G), then |v - r h v| m = 0(p 2 m ), 

0 < m < 2 

- for the quadratic elements (triangles, rectangles, tetrahedrons, 3-D 
rectangles) if: 

k = 2, v€H 3 (n) then, |v - r h v| m Q = 0(p 3 “), 

0 < m < 3 

- for the cubic elements (triangles, rectangles, tetrahedrons, 3-D 

rectangles) if: 

k = 3, v€H 4 (ft), then |v - r h v| m * 0(p* m ) 

0 < m < 4. 


Concerning the approximation error, they are optimal (i.e., the approxi- 
raation error is of the order of the best interpolation error), for instance, 
with the elements above, for Problem (3.1) - (3.5) when Yq = Y , = 0 

(Dirichlet problem) or Fq = 0, F^ » Y (Neuman problem), and ft is a 

polygon fully covered by the elements of the triangulation 
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4.b. Nonlinear Problems 

For nonlinear problems the situation is more difficult and the results 
are less complete. Usually convergence results can be proved by using energy 
type inequalities and convergence techniques which are appropriate for the 
type of equations considered: see for instance [28] for the nonlinear problem 
(3.13), (3.2), (3.3), and [29] for the Navier-Stokes equations. When 
compactness methods are used some involved compactness arguments for finite 
elements may be necessary: cf. in R. Temam [29] the proof of convergence of 
the nonconforming finite element methods for the Navier Stokes 
equations. Also by lack of uniqueness for nonlinear elliptic problems the 
convergence may be limited to a subsequence or may assume as usual that we are 
"close" to the solution. 

Error estimates are also more difficult to obtain than in the linear 
case. They assume usually more regularity on the equation and/or the solution 
that is necessary for convergence. 

5. Mixed and Hybrid Finite Elements 


5. a. Minimax Formulation of a Boundary Value Problem 

Consider an abstract boundary value problem of the form (3.12) 


( To find u€V such that 


a(u,v) = <£,v>, for all ueV 
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When the bilinear form a is futhermore symmetric, then (5.1) is equivalent 
to a convex minimization problem: 


(5.2) 


To minimize for v€V, 


J(v) = j a( v, v) - <£,v>. 


The infimum of J on V is attained at a unique point of V which is called 
a solution (or a minimizer) for the variational problem (5.2). In fact the 
solution of (5.1) is the same as that of (5.2). 

The mixed finite elements are closely related to duality. A natural 
framework for both questions arises when V is a linear subspace of a Hilbert 
space X of the form 


(5.3) 


V = {vex, b(v,<j> ) = 0 for all $€Y} , 


where Y is another Hilbert space and b is a bilinear continuous form on X 
x Y. We assume furthermore that a is extended as a bilinear continuous form 
on X and that Z is extended as a linear continuous form on X. 

In this case we introduce the Lagrangian of the problem (cf. Ekeland- 
Teraam [11]): 


(5.3) 


L(v,i|j) = J(v) + bCv, 1 ?). 


It is easily verified that 
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j J(v) if v€V 
Sup L(v ,<|> ) = \ 

<|>CY ( +» if v€X\V 

and that the minimization problem (for veV): 

(5.4) Inf { Sup L(v,*)} 

v£V t|»€ Y 

has the same solution and the same infimum as (5.2). 

Now we can associate with (5.3) the so-called dual problem of (5.4) which 
is a maximization problem in Y 


(5.5) Sup {inf L(v.’f)} 

ye Y v€V 

It is shown in [11] that if L (i.e., here b) satisfies a suitable condition, 
then (5.5) has a unique solution denoted <f> . Furthermore, the pair 
{u,(j>}ex x Y is a solution of (5.5) and (5.4) (or (5.1)) if and only if 
(u ,<j> ) = (u,4>) = 0, i.e.. 


(5.6) 


ja(u,v) + b(v ,<|> ) = <£,v>, for all v£X 
lb(u,i|i) = 0, for all t|>€Y 


The initial problem (5.1) (5.2) is written in X as a constrained 

minimization problem 


| To minimize J(v) for v€X, subject to the constraint 
^b(v,ijj) = 0 for all iJi€Y. 


(5.7) 



-29- 


The above framework associates to the initial problem (5.1) (5.2) (5.7) an 
element $ of X which is the Lagrange multiplier for the constrained 
optimization problem (5.7). 

The necessary condition on b which guarantees the existence of <|> , 
the so-called inf-sup condition, was introduced independently by Babuska [2] 
and Brezzi [5] and reads 



Equivalently (5.8) means that the linear operator B from X into Y' 
defined by 


(5.9) <Bv,i|>> = b(v,iji), for all v€x, for all iJi€Y' 


is an isomorphism from the orthogonal of V in X onto Y' or that the 
adjoint B' of B which maps X' into Y is an isomorphism from X onto 
the polar set V° of V 


V° = {0£X', <0 »v> = 0, for all v€V} . 


The reader is referred for more details to the article of Brezzi in this 
volume. Note that the form (5.6) of the problem can be studied independently 
of the corresponding Lagrangian and variational problems and is suitable for 
several types of generalizations: 
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- Given a linear continuous form x on Y * we can replace the second 
equation (5.6) by 

b(u,i|>) = for all ijiCY 

- More important, the form a may be nonlinear with respect to its first 

argument u, and this corresponds to considering nonlinear partial 

differential equations, in particular the Navier-Stokes equations (see below) 
or monotone operators (see [24]). 

See also in [11] a different point of view for duality which includes 
(5.6) as a particular case. 

Remark 5.1 . - Let us mention also here the penalization of (5.6) which 
leads to consideration of the following problems 

! To find u £ € X, <J> £ € Y such that 
a(u e ,v) + b(v,<t> £ ) = <£,v>, V vCX 
-ec(p e ,<|i) + b(u e ,<J>) = 0, for all iJi€Y 

where c is a bilinear continuous coercive symmetric form on Y and 
e > 0 is a fixed positive parameter which is intended to tend to 0. A 
solution u e ,<f, e of (5.10) exists for every e > 0, and u e ,4> e converges 

to the solution u, <J> of (5.6) when e+0; see M. Bercovier [3]. 
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5.b» Examples 
Stokes Equations 

The Stokes equations provide one of the most typical examples where the 
framework (5.6) is suitable. Stokes problem is the problem (1.7) - (1.10) 

when U = 0 and the nonlinear terra (u»V)u is dropped. In this case 

(ft C if 1 , n = 3 or more generally n ^ 3): 


X = Hj(ft) n = {v€L 2 (ft) n , |^_£L 2 (ft) n , Vi, v = 0 on r} 


V = {v€X, div v = 0} 


Y = {<|>eL 2 (ft), / 4»(x) dx = 0} 
ft 

n 3u 9v. 

a(u,v) = v 1 dx 


b( v ,<t> ) = -/ (div v) <|> dx. 
ft 


It can be shown that (5.6) is equivalent to the Stokes problem 


j 

( -vAu + grad p = 

f 

in 


(5.10) < 

\ 

div u = 

0 

in 

ft 

| 

L u = 

0 

on 

r. 


The operator div is a surjection from V onto Y (see for instance R. Temam 
[29]), and it follows immediately that the Babuska-Brezzi condition (5.8) is 


satisfied 
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As indicated above, we can set the Navier-Stokes equations in the frame- 
work (5.6) with a replaced by a nonlinear form 


a(u,v) = a Q (u,v) + a^Ujv) 


n 3u 9v 

a 0 (u > v) - " ? . I mr w: 
i,j=i n j j 


a (u,v) = l J u y-i v .dx. 
1 i , j=l a i 3 


In that case cf. [6]. 

Dirichlet Problem 

The framework (5.6) applies to (3.5) and provides the dual of this 
problem (see Ekeland-Temam [11]). For simplicity we restrict ourselves to the 
case where = F, = 0, i.e., we consider the Dirichlet problem in fi. 

We set 

X = hJ(«) x L 2 (fi) n , Y = L 2 (a) n 

v = l u = { u o» u i} ex » u i = Srad u Q } 
and for u = {uq,u^}, v = {vq,v^}€X and 'F€Y: 
a(u,v) = / Uq ,v^dx + / u^ • v^ dx 
<£,v> = / fv^dx 
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b(v,t|>) = / (Vj - grad v Q )»i|>dx. 

Si 

We identify (5.6) with (3.5). Condition (5.8) is trivially satisfied. The 
dual (5.5) reads 

1 2 

/ To maximize ~ f j I’M dx, 

(5.11) J 

( for 4i€L 2 (Cl) n , div i|> + f ■ 0 in Cl. 


Biharmonic Problem 


We consider the problem 


(5.12) A 2 u - f in Cl, 

(5.13) u - 0, Jii - 0 on r. 


It is set in the form (5.6) with 


V = Hq(0) = {v€L 2 (C1), 
on r} 


|-Xg— £L 2 (C1), for all i,j, and v = 0, = 0 

X i j 


a(u,v) = J Au Av dx 
Cl 


<Z,v> = / fvdx. 
Cl 


Then we set it in the form (5.6) with 
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X = x L 2 (Q), Y = L 2 (Q) 


V “ { u " { u o> u !} eX > u i = Au 0 l 


and 


for u = {u Q> u 1 } > v = { v q > v 1 }eX and i|»€Y 


a(u,v) = / u v dx, 

a 


<£,v> = J fv-dx 


b(v,ij>) = / (AVq - v ^ )»i|<dx. 


We identify (5.6) with (3.5). Condition (5.8) is trivially satisfied. The 
dual (5.5) reads 


(5.14) 


To maximize “ y / |*H^dx 


a 


for 4>£L (ft), Ai|i = f in ft. 


Problems involving the biharmonic appear in elasticity and in fluid mechanics 
for the treatment of the Stokes problem by utilization of a stream function. 


5.c. Mixed and Hybrid Elements 

Once we have reduced Problem (5.1) to (5.6), we are naturally led to 
approximate this last problem, i.e., 

- To find X^, Y^, a^, b^, c^ which approximate X, Y, a, b, c. 

- Solve for each h a discrete problem similar to (5.6): 
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To find u €X , <t>,£Y, such that 
h h h h 

VvV + VW " a h’V for a11 V 

b h^ u h’^h^ ~ for ^h eY * 

Finite element methods appear naturally in the construction of the 

spaces and Y^. We have more flexibility than in an ordinary finite 

element method since we can combine various finite elements for X, and 

h 

for Y^. The major difficulty arises in the verification of the condition 
(5.8) which leads sometimes to delicate algebraic questions. A thorough 
investigation of the inf sup condition for various finite elements related to 
the Navier-Stokes equations can be found in J. T. Oden and 0. P. Jacquotte 
[21]. In some cases, the number 8 in (5.8) corresponding to the discrete 
case depends on h and tends to 0 as h+0. In other cases, the inf sup 
condition does not hold in the discrete case and we can make it true by 
reducing the space Y^ in order to suppress the kernel of the discrete 
analogue of B^; in practice this amounts to a filtering procedure. 

The most famous example is the classical checkerboard instability for a Stokes 
problem corresponding to the Qj - Pg element: is a Q 1 approximation 

of H^(Q)^(n = 2) and Y^ is a Pg approximation of L^(ft); the 
filtering procedure is standard in this case; see also the analysis of the inf 
sup condition in Boland-Nicolaides [4] who show with a counter example that 
the best value of 8 in this (discrete) case is of the form c h. 

We will not develop further this question here since it is the object of 



[6] and other articles in this volume. 
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SPECIFIC APPLICATIONS 

Finite element methods have been the object of many applications in 
mechanics and physics. We describe now some specific applications. 

6. Navier-Stokes Equations 

The notations being the same as in Section 1 we consider the time 
dependent Navier-Stokes equations for viscous incompressible flow in a 
domain ft 

(6.1) - vAu + (u* V )u + Vp = f in ft x (0,T) 

(6.2) div u = 0 in ft x (0,T) 

(6.3) u — U on r x (0,T). 

The unknowns are the velocity vector u = u(x,t) and the pressure 
p = p(x,t); the voluraic forces f and the boundary velocity U (which may 
both depend on time) are given. 

We recall that from the mathematical point of view the initial value 
problem for the Navier-Stokes equations, i.e., (6.1) - (6.3) supplemented by 
an initial condition 

(6.4) u(x,0) = u Q (x), x£ft, u Q given, 

is well set in space dimension 2 (ft Cl?). However, we do not know yet if 
the same result is true in space dimension 3 (ft Cl?), i.e., we do not know if 
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the curl vector remains bounded or may become infinite even if the data are 
regular; see for instance R. Temam [29], 

The interval of time that we consider may be finite or infinite* Finite 
intervals of time occur naturally in the study of transient phenomena, while 
"infinite" intervals of time appear in the study of permanent regimes. For 
instance if f and U are independent of time, then in some cases the solu- 
tion u, p of (6.1) - (6.4) converges as t*« to a stationary solution, 
i.e., a solution of (1.7) - (1.9) . A sufficient condition for this to occur 
is that the Reynolds number R g is sufficiently small 


where U* is a typical velocity of the flow and L* a typical length of 
2 

ft • If R e is large, the convergence to a stationary solution is not 

guaranteed anymore. Based on experimental observations relative to turbu- 
lence, we actually expect that u(*,t), p(* ,t) do not converge anymore to 
time independent solutions even if the data f , U are independent of time . 

From the numerical point of view this will be the source of new difficulties 
which have not yet been explored and will not be addressed here (R. Temam 
[30]). Actually the computing power that is presently available leaves us at 
the threshold of the occurrence of nonstationary phenomena at least in space 
dimension two. 


9 Ct 0 1 

^Typical velocities are provided by V and f in the form L*lv 2 norra(U), 
norm(f), where appropriate norms are considered and the cl., 0 . are 
such that the corresponding expressions have the dimension of a velocity. 1 The 


sura of these two velocities is an appropriate definition of U*. 
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Up to now roost of the numerical computations on the full Navier-Stokes 
equations dealt with stationary phenomena (and transient phenomena). At this 
level a major difficulty for the numerical solution of (6.1) - (6.4) (or (1.7) 
- (1.9)) is the handling of the free divergence conditions which introduce 
complicated algebraic conditions in the discrete problem if it is not treated 
correctly. A certain number of methods related to or independent from the 
finite elements have been proposed to overcome the difficulties associated 
with the condition div u = 0. 

a) Utilization of the Penalty Method 
The penalty method which is due to R. Courant [10] in the context of con- 
strained optimization was applied to the Navier-Stokes equations in R. Temam 
[25], [26]. The idea which stems from the variational form of the Stokes 

problem (see Section 5) is to treat the condition div u = 0 as a constraint 
and to "penalize" it, i.e., to replace (6.1)(6.2) by 

3 u_ i 

(6.5) *-=- " vAu + (u »V)u - - V(V-u ) - f in fix (0,T) 

at e e s e e 

where e > 0 is a small parameter which is intended to tend to 0. It can 
be proved [26] that the solution of (6.3) - (6.5) converges to that of (6.1) - 
(6.4) when e+0. A full asymptotic expansion of u s , p e in terms of 
e can even be obtained in the simpler case of Stokes flows [29]. 

The penalty method has been applied in several ways by many authors to 
the finite element approximations of the Navier-Stokes equations, in 
particular with the mixed finite elements (see Remark 5.1). 
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b) Utilization of Algorithms: Artificial Time Dependence 

This method applies to the solution of the stationary problem (1.7) - 
(1.9) or to the solution of the stationary problems arising from time dis- 
cretization of (6.1) - (6.4). In these cases, an artificial evolution problem 
is introduced whose stationary solution is also the solution of the stationary 
Navier-Stok.es equation. 

For instance, one can consider the artificial evolution problem 


( 6 . 6 ) 


-vAu + (u*V)u + grad p = f 


+ a div u = 0 


(a > 0) or the equations of slightly compressible fluids 


(6.7) 


9 u 


- vAu + (u»V)u + grad p = f 


+ a div u = 0. 


In both cases the condition div u = 0 is not imposed at all times and 

follows simply from the properties of (6.6) and (6.7) for large t. 
Consequently, ordinary finite elements are used, i.e., finite elements not 
containing the condition div u = 0. 


c) Utilization of the Projection Method 
This method introduced in A. J. Chorin [8] and R. Teraara [27] is connected 
to the fractional step method. It consists in solving the time evolution of 
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(6,1) without (6.2) and then, more or less frequently, imposing (6.2) by 
projecting the velocity obtained on the free divergence vector fields. 

The time discretization (time mesh = At ) when U = 0 is given by 


( 6 . 8 ) 


/ 

m-l 

1 u 

- u 

J 

At 

1 

TT 111 

V u 

= u 

of 

u 


- vAu + (u *V)u - f in ft 


and u m = Proj. of u m which amounts to saying that 


m ~m , m 

u = u - grad q in ft 


(6.9) < div u m « 0 in 0 


V u m .v (= normal component of u m on T) m 0 on T 


Alternatively setting q m = At Ti m we can rewrite (6.9) as 


m 

u ~ u , m r\ 4 ^ 

— j-g + grad tt =0 in ft 


( 6 . 10 ) 


div u m = 0 in ft 


m ^ _ 

u = 0 on f , 


and p m provides an approximation for the pressure. It is, however, a poor 
approximation since it satisfies the following nonphysical boundary condition 
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that we infer from (6.10): 


( 6 . 11 ) 



on r . 


In order to determine u ra it is necessary to actually compute q m , ir ra , or 

at least their gradient; ir m is a solution of the Neuman problem which 
consists of (6.11) and 


( 6 . 12 ) 


A m 1 ~ra 

Air = — div u . 
At 


It seems better, for a more accurate determination of the pressure which 
avoids the undesirable boundary layer resulting from (6.11), to consider 
q m , rr ra as auxiliary functions and to compute the approximation p ra of the 
pressure by using the boundary value problem 


(6.13) p = Mu) P m = <Ku m ) 

that we deduce directly from the Navier-Stokes equations; cf. [29], At this 
point one can either use a Dirichlet or a Neuman boundary condition for p 
[15]. 

Many other forms of (6.8)(6.9) can be also considered: one can split the 
operators differently, leaving for example some viscosity in (6.9), one can 
use an explicit scheme in (6.8), or one can solve for the evolution (6.8) for 
several steps and perform the projection (6.9) periodically only. 

In all cases when the projection method is used we need a space of free 
divergence vector functions so that the projection (6.9) can be performed in a 
satisfactory manner. 
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d) Utilization of free divergence finite element spaces 

The simplest element, the piecewise linear (P^) function on triangles, 
cannot be used since the condition div u = 0 imposed on each triangle leads 
to too many algebraic relations and the spaces of discrete divergence free 
P^ vector functions may be reduced to the function 0, One can either impose 
the condition div u = 0 "less often," or go to more complicated elements 
such as the nonconforming P^ element (F. Thomasset) or P 2 , Qj, Q 2 > . .., 

elements. 

7. Fluid Structure Interactions 

In many industrial fields of interest, including the space industry, 
engineers are confronted with fluids (water, oil, kerosene, gases, ...) 
interacting with structures (tanks, containers, obstacles, ...) along a more 
or less extended area. 

In some cases deformations of the structure may be fairly important and 
even affect the motion of the fluid; engineers have then to solve problems 
including a coupling between fluid displacements and elastic deformations of 
the structure. 

We describe here the interaction between a free surface fluid and the 
structure, assumed to be elastic, which contains it, in an external force 
field. We follow J. Mathieu [13] and J. Mathieu, et al. [20], who computed 
the transient simulation of such a process with the BACCHUS Code. 

7,a. The Arbitrary Lagrange-Euler (ALE) Description 

■p 

We consider a moving domain ft (t) deforming with velocity 


w = w(x,t), and filled with an incompressible fluid of (constant) density 
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p 

p which obeys the Navier-Stokes equations. We consider also a second 

q F 

domain ft (t) made of elastic material and limiting ft . The fluid is 

g 

limited by a free surface S and the contact surface tt - ir(t) with ft* % 

The equations are 


P F {|_ v F + ((v F - w)V)v F } = div 0 F + p F f F in Q. 1 


div v =0 


(7.1) 


S 9 v .. S . S,S . n S 
p — = div ct + p f in 0 

O t 


^ (P S |J|) - 0 


where 


+ (w*V)p is the convection derivative associated with the 

fit 3 1 

vector field w 

v* = the velocity field in ft* 

a* = the stress tensor in ft 1 

i i Fx 

p = density in ft (constant in ft ) 

f « external forces 

S S 

J = Jacobian of the mapping ft (0)*ft (t). 
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The ALE description is determined by the actual velocity field 
is defined as follows: 

w = 0 in an internal eulerian region 

■n» 

w = v r on the free surface S(t) 

c 

w = v on the wet part of the wall ir(t). 



Figure 7 . 1 


The constitutive laws are 


o F = v [Vv F + (Vv F ) C ] + pi 

(W /V /s/ 


for the fluid, v = the kinematic viscosity, p = the hydrostatic 


a 


s (t) -*££) 

p S (o) 


F(t) T(t) YU)* 


w which 


pressure, 


and 
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for the solid where t Is the standard second Piola Kirchhoff tensor, F(t) 

A/ 

is the gradient of the mapping ^ ( 0 ^ ( t ) , and x(0) = a^(0) = 0. 

The boundary conditions are as follows : 

- On the structure , displacements are given on some part 

O 

T (t) of 3ft (t)Yit(t) and normal stresses are given along the remaining 

part of 3ft^(t)\ir (t) . 

- On the wet part of the wall , it ( t ) , we have no cavitation, 

F S 

V »v = v »v 

v the normal on 

along the wall 

and finally the normal stresses are continuous 

, F S. n 

(a - o )*v = 0. 

The discretization of the problem is made with a finite element 
discretization in space, providing an easy handling of the complex geometric 
configurations which are caused by large displacements. The elements are of 
degree one for the velocities (= for triangles, Qj for rectangles, 
transported for quadrilaterals), and piecewise constants for the hydro- 
static pressure. 

A one-step explicit finite difference scheme is used in time. This 
scheme is subject to the usual stability conditions limiting the time step. 


and we have a partial slip condition of the fluid 


°tg - - e(v t g - v t g > 
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The criteria taken into account are 

- free surface wave and viscous wave stability in the fluid 

- acoustic wave stability in the structure. 

Because of the more drastic limitation due to the stability criterion in the 

structure, a subcycling procedure is used for the structure calculation. 

A mesh adaptation procedure is necessary. At each fluid-calculation 

step, the free surface has to be repositioned. The displacements of the free 

surface and of the interface between the fluid and the structure induce a 

modification of the mesh in the mixed Lagrange-Euler region and possibly a 

modification of the Euler region (when its boundary intersects the free 

surface or the wall) or a degeneracy of some element. Thus a rezoning is 

automatically performed during the calculation. The rezoning is also 

necessary on ir(t) since any node belonging to 9ft (t)09ft (t) should be 

g 

at the same time a vertex of some element of ft (t) and of some element 
of ft F (t). 


A sample two-dimensional calculation is shown on Figure 7.2: 
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t = 0 D s and S are given. 

The internal mesh will be 
automatically designed. 



t = 0,2 The bottom of the tank is in a 
compressive phasis. 

Notice the behavior of the 
velocity fields on it: 
double-valued with normal 
component continuity. 



t g 0,3 


The bottom of the tank has 
entered an expansive phasis. 



t =* 0,5 The tank has reached large 
enough deformations. 

The free-surface tends to 
stabilize towards the 
horizontal. Notice the 
modification of the mesh in D^. 


Figure 7.2 
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8 . Thr ee Dimensional Euler Equations 

The finite element method has been applied to the computation of the 
solution of the Euler equations. We describe here the latest results in C. 
H. Bruneau, et al. [7], which concern the computation of steady vortex flows 
past a flat plate at high angle of attack. 

8.1. Description of the Problem 

The purpose of the computations in [7] is to investigate the developments 
of vortices at the tip of the plate and their propagation after the trailing 
edge. In the computations, the plate has no thickness and we expect a strong 
vortex structure to develop at the tip of the plate; it is not possible to use 
potential equations and the full Euler equations are necessary. For the 
computations, the plate is imbedded in a 3-D rectangular domain as shown on 
Figure 8.1; the aspect ratio is 0.5 and the angle of incidence is a = 15° 
or a = 30°. Figure 8.1 shows only half of the plate since symmetry with 
respect to the y variable is assumed. 

The incoming flow is given by q^ = (u^, v^ , w^) = (q coso, o, q 

i* 2 ^ 

2 

sina), a as above, and M the Mach number at infinity, M^ = — ^ , 

2 ° a 

o < 1« 

a = (y - 1)(H - . In these computations, the flow is subsonic, 

00 L 

M - 0.7 . 

00 
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Euler Equations 

They are written in conservative form 


3pu 3pv 3pw _ 

3x 3y TT~ U 


Conservation of mass 


( 8 . 1 ) 


3pu +p_ t 3puv^3puw 
<Fx 3y 3 z 


3puv 3pv + p 3pvw 
3 x <Ty 3 z 


= 0 


= 0 Conservation of momentum 


3puw , 3pvw , 3pw + p _ A 

~TT + ~TT + — Ti 0 


( 8 . 2 ) 


1 P + q _ H 

(y - Dp + T H 


Bernoulli's equation 


s 

where p, q = (u,v,w), p represents respectively 
velocity vector, and the pressure; H is the total and 
specific heats. 


the density, the 
Y the ratio of 


l 

I 
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Bouadary Conditions 

The boundary conditions on the plate are easy; the tangency condition on 
the plate (plane z = 0) reads w = 0, and the symmetry condition is v = 0 
on y = 0. On the contrary, the boundary conditions at the limits of the 
domain are not easy, especially downstream where it should allow the vortex to 
go through the exit plane. 

The flow variables have been set to the freestream values at the incoming 
boundaries (x = xq, z = zq), and no condition was imposed at the exit 
boundaries (x « Xj , z = z^) so that the values there are computed from the 
variational formulation. Due to the utilization of a least square method, 
this amounts simply to requiring that the first order equations be satisfied 
at those boundaries (see [7] for a discussion of the boundary conditions), 
y - yj the far field conditions v = 0 is used. 

8.2 The Numerical Method and the Numerical Results 

The equations are solved iteratively as follows: 

- A fixed point algorithm is based on equation (8.2), computing the 
density when u, v, w, p are known. 

- The nonlinear system (8.1) provides u, v, w, p, using the value of 

p from the previous step of the iteration. This system is linearized by 
Newton's method with linearized variables u, v, w, p. 

- The system for u, v, w, p is discretized by an appropriate 
finite element method (discretization of the conservative variables 
pu, pv, pw, ...) and the system is solved by a least square method. 

Figures 8.2 and 8.3 show a sample of the discretization grid. Finally, 
Figures 8.4 to 8.6 show the cross flow velocities at 70%, 90%, and 110% of the 


plate 






nnmmumui 
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FUTURE DEVELOPMENTS 

The advantages of FEM are well known. In particular 

- accurate and automating fitting of complicated geometries. 

- lowering (dividing by 2) the order of the differential operators that 
appear, due to the utilization of weak formulations. 

- avoiding the discretization of the boundary conditions of the Neuman 
type which disappear in the weak formulations. 

The inconveniences are 

- the need of a triangulation program for domains which has to be written 
once for all but is a discouraging preliminary step for the nonexpert. 

- the computational costs which by no way can compete with the achieve- 
ments of multigrid and spectral methods. 

However, the multigrid and spectral methods attain their optimal 
efficiency for rectangular domains. It is then conceivable that a combination 
of FEM and multigrid and spectral methods in relation with domain decomposi- 
tion and parallel computation can prove to be efficient. 

9. Domain Decomposition: Remarks on Future Developments 

The domain decomposition which lies at the foundation of the FEM can 
appear to be, in a different form, directly related to future developments in 
the method. Besides the geometrical considerations, there are many other good 
reasons for decomposing the solution of a boundary value problem in a large 
domain into the solution of similar problems in subdomains. Let us mention 


some of them: 
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a) Adaptive meshes 

Adaptive meshes may be suitable for geometrical or analytical reasons. 
For instance, the fitting of a curved domain with elements of different nature 
which are more appropriate in various parts of the domain. The refinement of 
the mesh in regions where the solutions are singular or nearly so (boundary 
layers, shocks, front flames, etc.) allow for extra computational effort to be 
concentrated in such regions. Thus finite elements of different types and 
sizes may be used in various parts of a domain, and a different treatment of 
the different regions can be useful. 

b) Physical Motivations 

Physical phenomena of a different nature may occur in different regions, 
and these regions should not be treated in a similar manner. For instance in 
aeronautics 

- a turbulence model is necessary near the airfoil 

- the Navier-Stokes equations with viscous effects and without turbulence 
are necessary at a certain distance but not too far 

-the Euler equations (i.e., no viscous effect) are sufficient far from the 
airfoil. 

Similar situations occur in combustion or in solid mechanics and lead 
naturally to the decomposition of the whole domain into smaller ones. 

c) Parallel Computation 

The technology of computers will apparently move more and more toward 
parallel computation. The decomposition of a domain into subdomains on which 
the computations are made simultaneously is appropriate for parallel computa- 
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tions. The difficulty is then to periodically synthesize the information from 
all the subdomains to ensure a correct interaction of the subdomains. 



Figure 9.1 
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technique is satisfactorily mastered, then one can consider, as has already 
been done (see for instance [ 13] [23] [ 19] ) , combining the advantages of FEM and 
spectral and multigrid methods by using domain decomposition as suggested in 
Figure 9.1. This has already been done, but the goal is to get as close as 
possible to the performance of the fast methods. 
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CONCLUDING REMARK 

The FEM cannot pretend to be the best method in all situations, but it 
has proved to be an efficient and performant method in many cases and can hope 
to be the object of future interesting developments. R. Feyman says in his 
book [12] that he was able to solve some problems that other people could not 
solve just because he had some tools in his box that others did not have. It 
would be too bad not to have the FEM tool in his box of numerical methods. 


-57- 


REFERENCES 


[1] J. H. Argyris, "Energy theorems and structural analysis, part I: 

General Theory," Aircraft Engineering , 26, 1954, pp. 347-356, 383-387, 
394, 27, 1955, pp. 42-58, 80-94, 125-134. 

[2] I. Babuska, "The finite element method with Lagrangian multipliers," 

Numer. Math. , 20, 1973, pp. 179-192. 

[3] M. Bercovier, "Perturbation of mixed variational problems," Application 
to mixed finite element methods, RAIRO Anal. Num. , 12, 1978, pp. 211— 
236. 

[4] J. Boland and R. Nicolaides, "Stable and semistable low order finite 

elements for viscous flows," University of Connecticut, preprint, 1984. 

[5] F. Brezzi , "On the existence, uniqueness and approximation of saddle- 

point problems arising from Lagrangian Multipliers," RAIRO-Analyse 

Numerique , 8, 1974, pp. 129-151. 

[6] F. Brezzi, "Finite element approximations of the Von KarraSn equations," 
RAIRO Numer. Anal., 12, No. 4, 1978, pp. 303-312. 



-58- 


[7] C. H. Bruneau, J. J. Chattot, J. Larainie, and R. Temam, "Computation of 
vortex flows past a flat plate at high angle of attack," in Proceedings 
of the 10th International Conference on Numerical Methods in Fluid 
Mechanics (Peking 1986), Lecture Notes in Physics, Springer-Verlag, New 
York, 1987. 

[8] A. J. Chorin, "A numerical method for solving incompressible viscous 

flow problems," J . Comp . Phys . , 2, 1967, pp. 12-26. 

[9] P. G. Ciarlet, The Finite Element Method for Elliptic Problems , North- 

Holland, Amsterdam, 1978. 

[10] R. Courant, "Variational methods for the solution of problems of 

equilibrium and vibrations," Bull. Amer. Math. Soc. , 49, 1943, pp. 1-23. 

[11] I. Ekeland and R. Temam, Convex Analysis and Variational Problems , 

Nor th-Hol land, Amsterdam, 1976. 

[12] R. Feyman, Surely You"re Joking Mr. Feyman , W. W. Norton & Company, New 
York, 1985. 

[13] B. Gest, "Methode de resolution rapide des systemes lin€aires," Thesis, 
Laboratoire d'Analyse Num£rique University Paris Sud, Orsay, 1985. 

[14] V. Girault and P. A. Raviart, Finite element methods for the Navier- 
Stokes equations, Springer Series in Computational Mathematics, 1986. 


-59- 


[15] P. M. Gresho and R. L. Sani, "On pressure boundary conditions for the 
incompressible Navier-Stokes equations," University of Colorado, 
Boulder, preprint, 1986. 

[16] J. Leray, "Etude de diverses Equations integrales non lineaires et de 
quelques problemes que pose l'hydrodynamique," J. Math. Pures Appl. , 12, 
1933, pp. 1-82. 

[17] J. Leray, "Essais sur les mouvements plans d'un liquide visqueux que 
limite des parois," J. Math. Pures Appl. , 13, 1934, pp. 331-418. 

[18] J. Leray, "Essai sur le mouvement d'un liquide emplissant l'espace," 
Acta Math ., 53 , 1934, pp. 193-248. 

[19] C. Lesartre, Thesis in preparation, Laboratoire d' Analyse Num&rique, 
Universite' Paris Sud, Orsay, 1987. 

[20] J. Mathieu, P. Ravier, J. Boujot, P. Gendre, M. Hittinger, "Interaction 
between structure and free surface fluid with large displacements of by 
finite elements," in Proceedings of the 10th International Conference on 
Numerical Methods in Fluid Mechanics (Peking 1986), Lecture Notes in 
Physics, Springer-Verlag, New York, 1987. 

[21] J. T. Oden and 0. P. Jacquotte, "Stability of some mixed finite element 
methods for Stokesian flows," Computer Methods in Applied Mechanics and 
Engineering, 43 , 1984, pp. 231-248. 



-60- 


[22] J. Periaux, First International Conference on Industrial and Applied 
Mathematics, Paris, June 1987. 

[23] B. Ravier, "Methode de decomposition des domaine," Thesis, Laboratoire 
d'Analyse Numerique, University Paris Sud, Orsay, 1986. 

[24] B. Scheurer, "Existence et Approximation de points selles pour certains 
problemes nonlin§aires," RAIRO-Analyse Numerique , 11 , 1977, pp. 369—400. 

[25] R. Temam, "Sur l'approximation des solutions des Equations de Navier- 
Stokes," C. R. Ac. Sc. Paris , S€rie A, 262, 1966, pp. 219-221. 

[26] R. Temam, "Une inethode d'approximation de la solution des Equations de 
Navier-Stokes ," Bull. Soc. Math. France , 98 , 1968, pp. 115-152. 

[27] R. Temam, "Sur l'approximation de la solution des Equations de Navier- 
Stokes par la methode des pas f ractionnaires II," Arch. Rat. Mech. 
Anal. , 33 , 1969, pp. 377-385. 

[28] R. Temam, Numerical Analysis , Reidel Pub. Comp., Dordrecht, Holland, 
1973. 

[29] R. Temam, Navier-Stokes Equations , 3rd revised edition, North-Holland, 


Amsterdam, 1984 


- 61 - 


[30] R. Temam, "Chaos in Fluids," Lecture at the SIAM National Meeting, 
Boston, July 1986. 

[31] R. Teraam, Lecture at the INRIA-SIAM-S0CA1 Workshop on Domain 
Decomposition, Paris, Proceedings to be published by SIAM, Philadelphia, 
January 1987. 

[32] F. Thomasset, Implementation of Finite Element Methods for Navier-Stokes 
Equations , Springer Series in Comp. Phys., 1981. 

[33] 0. C. Zienckiewicz, The Finite Element Method in Engineering Science , 
MacGtaw-Hill, London, 1971. 



Standard Bibliographic Page 


1. Report No. NASA CR- 178222 2. Government Accession No. 

ICASE Report No. 86-77 

3. Recipient’s Catalog No. 

4. Title and Subtitle 

SURVEY OF THE STATUS OF FINITE ELEMENT METHODS 
FOR PARTIAL DIFFERENTIAL EQUATIONS 

5. Report Date 

November 1986 

6. Performing Organization Code 

7. Author (s) 

Roger Teraara 

8. Performing Organization Report No. 

86-77 

9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665-5225 

10. Work Unit No. 

11. Contract or Grant No. 

NAS1-18107 

13. Type of Report and Period Covered 

Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

sns-Qn-9i-ni 

15. Supplementary Notes 

Langley Technical Monitor: 
J. C. South 

Final Report 


16. Abstract 

The finite element methods (FEM) have proved to be a powerful technique for 
the solution of boundary value problems associated with partial differential 
equations of either elliptic, parabolic, or hyperbolic type. They also have a 
good potential for utilization on parallel computers particularly in relation to 
the concept of domain decomposition. 

This report is intended as an introduction to the FEM for the nonspecialist. 

I It contains a survey which is totally nonexhaustive, and it also contains as an 
j illustration, a report on some new results concerning two specific applications, 
namely a free boundary fluid-structure interaction problem and the Euler equations 
for inviscid flows. 


17. Key Words (Suggested by Authors (s)) 

finite elements, Euler equations 
structural mechanics 


18. Distribution Statement 

64 - Numerical Analysis 
05 - Aircraft Design, Testing, 
and Performance 

Unclassified - unlimited 


19. Security Classif.(of this report) 

20. Security Classif.(of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 

63 

A04 


For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley, 1987 


